Equilibration through local information exchange in networks 
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We study the equilibrium states of energy functions involving a large set of real variables, defined 
on the links of sparsely connected networks, and interacting at the network nodes, using the cavity 
and replica methods. When applied to the representative problem of network resource allocation, 
an efficient distributed algorithm is devised, with simulations showing full agreement with theory. 
Scaling properties with the network connectivity and the resource availability are found. 
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Many theoretically challenging and practically impor- 
tant problems involve interacting variables connected by 
network structures [lj. Statistical mechanics of disor- 
dered systems makes contributions towards the under- 
standing of such systems at two levels. Macroscopically, 
it describes the typical behavior of the networks, using 
techniques such as the replica method. Microscopically, it 
analyzes the relation between variables, using techniques 
such as the cavity method, that give rise to efficient com- 
putational algorithms. In computer science, probabilistic 
inference based on graphical structures has been devel- 
oped and applied 00 , mainly for providing approximate 
solutions to specific instances of limited size. Examples of 
recent success included the belief propagation algorithm 
for error-correcting codes £| and the survey propagation 
algorithm for the satisfiability problem 0]. 

Most analyses so far have focused on networks of dis- 
crete variables. However, many typical problems, such as 
network resource allocation, involve continuous variables. 
Compared with discrete variables, analyses for continu- 
ous variables were much less explored. The main obstacle 
comes from the need to pass among the nodes entire free 
energy functions as messages. This is much more com- 
plex than cases of discrete values, where the messages 
are countable sets of conditional probability estimates 
of discrete values. Previous work in the computer sci- 
ence literature focused on modeling these functions for 
getting good approximations in feasible time scales 0. 
There have been attempts to simplify the messages for 
continuous variables, for example, to parametrize them 
using eigenfunction decomposition for special cases, but 
the general feasibility remains an open question 0. 

In this Rapid Communication we study a system of 
real variables that can be mapped onto a sparse graph. 
Based on the analysis, we demonstrate the close relation- 
ship between belief propagation algorithms and the Bethe 
approximation in statistical physics 0, and propose a 
message-passing approximation method, generally appli- 
cable to problems of continuous variables. The method 
is efficient since the messages consist of only the first and 
second derivatives of the vertex free energies derived from 
our analysis. The key to the successful simplification, not 
needed for the simpler case of discrete variables, is that 
the messages passed to a target node are accompanied by 



information-provision messages from the target node, to 
first determine the state at which the derivatives should 
be calculated. 

We first formulate the problem at a general temper- 
ature, and then focus on a prototype for optimization. 
The traditional approach for optimization on networks is 
to adopt computationally demanding global optimization 
techniques, such as linear or quadratic programming 0. 
In contrast, message-passing approaches have the poten- 
tial to solve global optimization problems via local up- 
dates, thereby reducing the computational complexity. 
An even more important advantage, relevant to practical 
implementation, is its distributive nature. Since it does 
not require a global optimizer, it is particularly suitable 
for distributive control in large or evolving networks. 

We consider a sparse network with N nodes, labeled 
i = l,...,N. Each node i is randomly connected to 
c other nodes. The connectivity matrix is given by 



Ai 



1,0 for connected and unconnected node pairs 



respectively. A link variable is defined on each con- 
nected link from j to i. We consider an energy function 
(cost) E = J2 {lj) ■■ljO{!Ji J )~Y,, ip(Xi, {Vij\Aij = 1}), where 
Xi is a quenched variable defined on node i. In the context 
of probabilistic inference, y.ij may represent the coupling 
between observables in nodes j and i, 4>{y.ij) may cor- 
respond to the logarithm of the prior distribution of yij , 
and ip(Xi, {i/ij \Aij — 1}) the logarithm of the likelihood of 
the observables A^. In the context of resource allocation, 
Uij = —yji may represent the current from node j to i, 
4>(yij) the transportation cost, and ip(Xi, {yij\Aij = 1}) 
the performance cost of the allocation task on node i, 
dependent on the node capacity A^. 

We are interested in the case of intensive connectivity 
c ~ 0(1) <c N. Since the probability of finding a loop 
of finite length on the network is low, the cavity method 
well describes the local environment of a node. A node 
is connected to c branches in a tree structure, and the 
correlations among the branches of the tree are neglected. 
In each branch, nodes are arranged in generations. A 
node is connected to an ancestor node of the previous 
generation, and other c— 1 descendent nodes of the next 
generation. Considering node i as the ancestor of node 
j, the descendents of node j form a tree structure T 
with c — 1 branches, labeled by k ^ i for Ajk = 1. At a 
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temperature T = ft ' , the free energy F(yij\T) can be 
expressed in terms of the free energies F(yj k \T k ) of its 
descendents. The free energy can be considered as the 
sum of two parts, F(y\T) = N T F av +F v (y\T), where N T 
is the number of nodes in the tree T, F av is the average 
free energy per node, and Fv(y\T) is referred to as the 
vertex free energy. This leads to the recursion relation 



-F av , (1) 



F V ( yij \T) = - rin |n (/ ex p 



p — 



A jk =l 



TUnj H (J dy jk j exp 



-Ptp{Xj,{yj k }) 



-pJ2( F v(yjk\T k )+cf>(y jk )) 



,(2) 



A ik =H A 



where T k is the tree terminated at node fc, and (•■■}> 
represents the average over the distribution of A. Inter- 
estingly, the recursive relation of Eq. Q can be directly 
linked to probabilistic message passing (belief propaga- 
tion), where the logarithms of messages passed between 
nodes are proportional to the vertex free energies. 

For more concrete discussions, we focus on a proto- 
type for optimization, termed resource allocation and 
well known in the areas of computer science and oper- 
ations management [lJJ, EH • The analysis of the problem 
is applicable to typical situations where a large number 
of nodes are required to balance loads and/or resources, 
such as reducing internet traffic congestion and stream- 
lining network flows of commodities [Tsi ]. In computer 
science, many practical solutions are usually heuristic 
and focus on practical aspects (e.g., communication pro- 
tocols) . Here we study a more generic version of the prob- 
lem represented by nodes of some computational power 
that should carry out tasks. Both computational powers 
and tasks will be chosen at random from some arbitrary 
distribution. The nodes are located on a randomly cho- 
sen sparse network of some connectivity. The goal is to 
allocate tasks on the network such that demands will be 
satisfied while the migration of (sub-)tasks is minimized. 

We focus here on the satisfiable case where the to- 
tal computing power is greater than the demand, and 
where the number of nodes involved is very large; the 
unsatisfiable case can be investigated using a similar ap- 
proach Each node on the network has a capacity 
(computational capability minus allocated tasks) ran- 
domly drawn from a distribution p(Aj). With the aim to 
satisfy the capacity constraints, we have if>(\i, {yij\Aij = 
l})=ln[6(- J2j Aj2/y-A;)+e], where e -> 0. The prob- 
lem reduces to the load balancing task of minimizing the 
energy function (cost) E — J2^j) ^ij 4>(Uij) ) subject to the 
capacity constraints ^ . Ai 3 yi 3 ; + K > 0. 



When <j){y) is a general even function of the current y, 
we may also derive Eq. using the replica method. We 
first introduce the chemical potentials fa of nodes i, and 
approximate the current y+j as driven by the potential 
differences between nodes yij = [ij—fii. Since sparse net- 
works are locally treelike, the probability of finding short 
loops is vanishing in large networks, and the approxima- 
tion works well. 

Considering the optimization problem in the space of 
chemical potentials, we calculate the replicated partition 
function {Z n )^\ averaged over the connectivity matrix 
and capacity distribution, and take the limit n — > 0. 
Assuming replica symmetry, the saddle point equations 
yield a recursion relation for a two-component function 
R dependent on the tree structure T, given by 

^,M|T) = ^n(/ d/i feJ R(/i,/i fe |T fe )] e(j2fx k 
k=i ^ J ' \k=\ 



-c:{i+z + \ v(T) j exp|-y^ 2 - (3 ^0" _ Mfc)J) ( 3 ) 

where I? is a constant, represents the tree terminated 
at the kth descendent, and Xy(T) the capacity of the 
vertex of the tree T. The term /3e/it 2 /2, with e — > 0, 
is introduced to break the translational symmetry of the 
chemical potentials, since the energy function is invariant 
under the addition of a constant to all chemical poten- 
tials. 

Equation J3J) expresses R(z, /i|T) in terms of c— 1 func- 
tions R((i, fi k \T k ) (k = l,..,c— 1), a characteristic of 
the tree structure. Furthermore, except for the factor 
exp(— /3e/x 2 /2), R is a function of y = fx — z, which is 
interpreted as the current drawn from a node with chem- 
ical potential by its ancestor with chemical potential z. 
One can then express the function R as the product of 
a vertex partition function Zy and a normalization fac- 
tor W, that is, R{z,fi\T) = W(z)Z v (y\T). In the limit 
e^0, the dependence on \i and y decouples, enabling one 
to derive a recursion relation for the vertex free energy 
F v (y\T) = -T\nZ v (y\T) and arrive at Eq. ©. 

The current distribution and the average free en- 
ergy per link can be derived by integrating the cur- 
rent y 1 in a link from one vertex to another, fed by 
the trees Ti and T 2 , respectively; the obtained expres- 
sions are P(y) = (5(y — y'))+ and (E) = {^(y'))* where 
(•)* - (Jdy'e X p[-f3E(y')}(*)/Jdy'exp[-(3E(y')}) x , 
and E(y')=F v (y>\T 1 )+F v (-y>\T 2 )+<f>(y'). 

Figure 1(a) shows results of the iteration of Eq. 
in the case of optimization (T = 0) based on discretizing 
Fy(y\T) as a vector. The capacity distribution p(X) is 
Gaussian of variance 1 and average (A). Each iteration 
corresponds to adding one extra generation (1000 new 
nodes in our simulations) to the tree structure, such that 
the iterative process corresponds to approximating the 
network by an increasingly extensive tree. We observe 
that after an initial rise with iteration steps, the average 
energies converge to steady-state values, at a rate which 
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FIG. 1: Results for TV = 1000 and cj>(y) = y' 2 /2. (a) (tj>) 
obtained by iterating Eq. Q as a function of t for (A) = 0.1, 
0.2, 0.4, 0.6, 0.8 (top to bottom) and c=3. Dashed line: The 
asymptotic {</)} for (A) =0.1. Inset: 7 as a function of (A), 
(b) K 2 ((f>) as a function of (A) for c = 3 (O), 4 (□), 5 (0), 
large K (line). Inset: K 2 (cf>) as a function of time for random 
sequential update of Eqs. J1J and 0. Symbols: same as (a). 
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FIG. 2: Results for N= 1000 and (f){y)=y 2 /2. (a) The current 
distribution P(Ky)/K for (A) =0.02,0.5, 1, and c = 3 (solid 
lines), 4 (dotted), 5 (dot-dashed), large K (long dashed). In- 
set: P(y = 0) as a function of (A) for c = 3 (O), 4 (□), 
5 (0), large K (line), (b) The resource distribution P(r) 
for (A) = 0.02, 0.1, 0.5, large K. Symbols: as in (a). Inset: 
P(r>0) as a function of (A). Symbols: as in the inset of (a). 



increases with the average capacity. 

To study the convergence rate of the iterations, we 
fit the average energy at iteration step t using (E(t) — 
E(oo)} ~exp(— jt) in the asymptotic regime. As shown 
in the inset of Fig. 1(a), the relaxation rate 7 increases 
with the average capacity. A cusp appears at the aver- 
age capacity of about 0.45, below which convergence is 
slow due to a plateau that develops in the average energy 
curve before the final stage. The slowdown is probably 
due to the appearance of increasingly large clusters of 
nodes with negative resources, which draw currents from 
increasingly extensive regions of nodes with excess re- 
sources to satisfy the demand. 

The local nature of the recursion relation points 
to the possibility that the network optimization can be 
solved by message-passing approaches. Instead of pass- 
ing the functions Fy(y\T) of the current y as messages, 
we simplify each message to two parameters, namely, the 
first and second derivatives of the vertex free energies. 
Let {AijiBu) = [dF v {y ij \T j )/dy ij ,d 2 F v {y ij \T j )/dyU 
be the message passed from node j to i. Using Eq. 0J, 
the recursion relations lead to the message (Aij , Bij ) 



Bi 



Q(-fiij + e) 



(4) 



Mij : 



+Xj - Vij 



A Jk[yjk - {<t>' jk + A jk ){4>% + B jk ) : ] 



,0 , 



(5) 



with 1 



and 1 



ujf. auu ijjjfc representing the first and second deriva- 
tives of 4>{y) at y = yj k , respectively. 

The algorithm is complete with the determination of 
the drawn current y^ at which the derivatives comprising 
the messages should be computed. Two methods are pro- 
posed. In the first, when messages are sent from node j 
to the ancestor node i, backward messages yj k computed 
from the same optimization step are sent from node j 
to the descendent nodes fc, informing them of the par- 
ticular arguments to be used for calculating subsequent 
messages. In the second, node j first receives the mes- 
sages (Aji, Bji) and current yji from the ancestor node i, 
and update the current by minimizing the total cost. 
Both methods work well for the quadratic cost functions. 

For comparison, an independent exact optimization is 
available at zero temperature. The chemical potentials 
turn out to be the Lagrange multipliers of the capacity 
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constraints, and the relation between the currents and 
the chemical potentials turns out to be exact. The Kiihn- 
Tucker conditions for the optimal solution yield 



1 



^Ai^j+Xi ,0 



(6) 



Like in the message-passing algorithm, this condition also 
provides a local iterative solution to the optimization 
problem. Simulations show that it yields excellent agree- 
ment with Eqs. CEJ, Q, and (J5J. 

To study the dependence on the connectivity, we 
first consider the limit of large K = c — 1 , in which 
Eq. (0} converges to the steady-state results of Ay = 
maxfJf- 1 J2k& -AjkAjk - 0] and B iS ~ K' 1 . Then 
J2k^i AjkAjk becomes self-averaging and equal to KrriA, 
where bia^-K" -1 is the mean of the messages Ay. Thus, 
Vij ~ Mi ~ K . The physical picture of this scaling 
behavior is that the current drawn by a node is shared 
among the K descendent nodes. After rescaling, quan- 
tities such as K 2 (</>), P(Ky)/K and P(Kfi)/K become 
purely dependent on the average capacity (A). 

For increasing finite values of K, Fig. 1(b) shows the 
common trend of K 2 (ip) decreasing with (A) exponen- 
tially, and gradually approaching the large K limit. The 
scaling property extends to the optimization dynamics 
[Fig. 1(b) inset]. As shown in Fig. 2(a), the current dis- 
tribution P(Ky)/K consists of a 5 function component 
at y = and a continuous component, whose breadth de- 
creases with (A). Remarkably, the distributions for differ- 
ent connectivities collapse almost perfectly after the cur- 
rents are rescaled by K , with a very mild dependence 
on K and gradually approaching the large K limit. As 
shown in the inset of Fig. 2(a), the fraction of idle links 
increases with (A). The fraction has a weak dependence 
on the connectivity, confirming the almost universal dis- 
tributions rescaled for different K . 

Since the current on a link scales as if -1 , the allo- 
cated resource of a node should have a weak dependence 
on the connectivity. Defining the resource at node i by 
n = Xi+^2j •Aijyij, the resource distribution P{r) shown 



in Fig. 2(b) confirms this behavior even at low connec- 
tivities. The fraction of nodes with unsaturated capacity 
constraints increases with the average capacity, and is 
weakly dependent on the connectivity [Fig. 2(b) inset]. 
Hence the saturated nodes form a percolating cluster at 
a low average capacity, and breaks into isolated clusters 
at a high average capacity. It is interesting to note that 
at the average capacity of 0.45, below which a plateau 
starts to develop in the relaxation rate of the recursion 
relation, Eq. QJ, the fraction of saturated nodes is about 
0.47, close to the percolation threshold of 0.5 for c=3. 

In summary, using the example of the resource alloca- 
tion problem on sparsely connected networks, we stud- 
ied the use of message-passing methods for equilibration 
using both replica and cavity based analyses. A local 
algorithm was devised and successfully applied to this 
task. The study also reveals the scaling properties of 
this model, showing that the resource distribution on the 
nodes depends principally on the networkwide availabil- 
ity of resources, and depends only weakly on the connec- 
tivity. Links share the task of resource provision, leading 
to current distributions that are almost universally de- 
pendent on the resource availability after rescaling. 

While the analysis focused on fixed connectivity and 
zero temperature, it can accommodate any connectivity 
profile and temperature parameter and may be used for 
analyzing a range of inference problems. For instance, 
we have considered the effects of adding anharmonic and 
frictional terms to the quadratic cost function. The 
message-passing function can be adapted to these vari- 
ations, and the results will be presented elsewhere |13[ . 
Both analysis and algorithm extend the use of current 
message-passing techniques to inference in problems with 
continuous variables, opening up a rich area for further 
investigations with many potential applications. 
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